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ABSTRACT 



The HD 196885 system is composed of a binary star and a planet orbiting the primary. The orbit of the binary is fully constrained 
by astrometry, but for the planet the inclination with respect to the plane of the sky and the longitude of the node are unknown. Here 
we perform a full analysis of the HD 196885 system by exploring the two free parameters of the planet and choosing different sets of 
angular variables. We find that the most likely configurations for the planet is either nearly coplanar orbits (prograde and retrograde), 
or highly inclined orbits near the Lidov-Kozai equilibrium points, i = 44° or i = 137°. Among coplanar orbits, the retrograde ones 
appear to be less chaotic, while for the orbits near the Lidov-Kozai equilibria, those around cj = 270° are more reliable, where a) is 
the argument of pericenter of the planet's orbit with respect to the binary's orbit. From the observer's point of view (plane of the sky) 
stable areas are restricted to (I\,0.\) ~ (65°, 80°), (65°, 260°), (1 15°, 80°), and (1 15°, 260°), where I\ is the inclination of the planet 
and f2i is the longitude of ascending node. 



Key words, celestial mechanics - planets and satellites: dynamical evolution and stability 
close - (Stars:) individual: HD 196885 



(Stars:) binaries (including multiple): 



1. Introduction 

About 20% of all known exoplan ets have b een fou nd to 
inhabit multiple stellar systems (Desidera & Barbieri 2007; 
Mugr auer & Neuhau ser 2009). The theories of planet forma- 
tion around a star that is part of a binary system have con- 
siderable challenges. A brief review of theories can be seen in 
lOuintana et"ai1 (120021) : iThebauld d201ll) . However, it is very im- 
portant to know the real masses and the spatial configuration of 
such systems in order to better understand the processes involved 
in their formation. 

Among planets in binaries, a few of them are found in com- 
pact systems with semimajor axis o f the binary les s than 1 00 AU: 
HD 196885 (Correia et al. 2008: IChauyin et alj l201lb. Gl 86 



(Oueloz et al. 2000; Lagrange et al. 2006), y-Cep (Hatzes et al 



2003 



2004) 



Neuhause r etalJ 120071) . "and HD41004 (Zuck er et al. 

Even for these tight systems, the poor precision of first 
epoch observations together with incomplete time span of the 
observations lead t o severa l best fi t solutions wi th almost the 
same residuals (e.g iTorresI (120071) : ICorreia et all d2008)). It is 
then advisable to keep in mind this wide range of possible con- 
figurations, and to combine radial velocity fits with pos sible for- 
mation scenarios (e.g. y-Cephei, Giupp one et al.1 1201 1). 

Recent observations of the HD 196885 system allowed to 
constrain the orientation of the binar y orbit combining astromet- 
ric and radial velocity observations (IChauvin et al.ll201lt) . Thus, 
the real mass of the stellar companion was esta blished at 0.45 
M(7). A few numerical integrations were carried bv lChauvin et al.l 
(1201 ll) . however it is not clear which are the compatible regions 
of solutions for the planet around the central star and why some 
individual coplanar solutions are unstable. We intend to clarify 
this picture by performing massive numerical integrations over 
the entire space of the free parameters and compare the results 
with analytical models. This will allow us to completely clar- 



ify all the possible dynamical regimes in the HD 196885 system, 
and to put constraints on the forthcoming observations. 



2. Best fit solution and reference angles 

HD 196885 B, the stellar companion of HD 196885 A, was de- 
tected combining imaging and spectroscopic observations, as be- 
ing a M1±IV dwarf located at 0.7", which corresponds to 23 AU 
in projected physical separation (Chauvin et al. 2007). 

The planetary companion, HD 196885 Ab, was detected 
using radial velocity data from ELODIE, C ORALIE and 
CORAVEL observations spread over 14 years (ICorreia et al.l 
2008). The orbital solution for the planet gave a minimum mass 
of niAb sini = 2.96 Mj up , a period of P = 3.69 + 0.03 yr, and 
an eccentricity of e — 0.462 + 0.026. The authors additionally 
constrained the binary companion HD 196885 B with a period 
of 40 < P < 120 yr, a semi-major axis 14 < a < 30 AU, 
and a minimum m ass of 0.3 < niR sin i < 0.6 M Q . Based on 
Lick observations, Fisc her et al.1 (12009) confirmed this solution 
for the system. In both studies, a relatively large range of masses 
and periods were found compatible with residuals for the stellar 
companion. 

More recently, Chauvin et al (201 1) combined five astromet- 
ric measurements obtained with VLT/NACO spread over 4 years 
and all sources of radial velocity data (four sets, summarizing 
187 observations) to determine the full orbit of the binary in the 
space and a consistent projection of planetary parameters. We 
show the best fit obtained by Chauvin et al. (201 1) in Table Q] 

The characterization of the orbits in the HD 196885 system 
by Chauvin et al. (201 1) left only two parameters undetermined, 
the inclination of the planet with respect to the plane of the sky, 
I\, and the longitude of the node, Qi (TabUJ. It is then possible to 
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Table 1. Orbital parameters for the HD 196885 system obtained 
by Chauvin et al. (201 1) for JD = 2455198. 



Param. 


[unit] 


orbit 1 (planet) 


orbit 2 (binary) 


di 


[AU] 


2.6 + 0.1 


21.00 ±0.86 


e, 




0.48 ± 0.02 


0.42 ± 0.03 


(i); 


[deg] 


93.2 ±3.0 


241.9 ±3.1 


M t 


[deg] 


349.1 ±1.80 


121 ±45 


Q., 


[deg] 


? 


79.8 ±0.1 


1, 


[deg] 


? 


116.8 ±0.7 


IUj 


[M Jup ] 


2.98/ sin /i 


472 




~<5rbitl 



Fig. 1. Fundamental planes for the definition of the reference an- 
gles. The angles in black are determined directly from the ob- 
servations (Tab.[TJi, while the two angles in blue (I\ , AQ) are the 
missing angles in the observations that prevent the full character- 
ization of the orbits in the system. The angles in red (i, Q, w) are 
also undetermined, but related with the other angles. They are 
linked to the physical system, thus independent of the reference 
frame, and more suitable to study the dynamics. 



(e.g. [Smart 1965). If one studies the system for mutually inclined 
configurations (/, Q), then we get for the remaining variables: 



cos I\ - cos I2 cos i - sin h sin i cos Q , 

(cos I2 - cos I\ cos i) 

cos Aw = , 

(sin I\ sin 1) 

cos AQ = cos Q cos Aw + sin Q sin Aw cos i , 



(1) 
(2) 
(3) 



noting that when 0° < Q < 180°, we have 0° < Aw < 180° 
and 0° < AQ < 180°, and for 180° < Q < 360°, we have 
180° < Aw < 360° and 180° < AQ < 360°. 

Alternatively, if one prefers to explore the space of possi- 
ble solutions starting from the observer's point of view (7i, Qj), 
the spherical triangle defined by the sides AQ, Aw, Q has three 
known parameters AQ, I\ , li and the following relations can be 
used to determine the physical system: 



cos; 



cos Aw 



cos I2 cos I\ + sin h sin 7) cos AQ . 
(cos I2 - cos I\ cos i) 



(sin/i sin/) 
cosQ = cos AQcos Aw - sinAQsinAwcos/i . 

3. Analytical model 



(4) 
(5) 
(6) 



Before studying the full massive problem, it is useful to look at 
the restricted inner problem (the orbit of the outer companion 
is fixed). This approximation is not very far from being true, 
since the mass of the planet is much smaller than the mass of 
the stellar bodies (m\ <k mo, nv£), and can therefore be seen as 
a "test particle" (;«i = 0). The restricted problem is easier to be 
studied and allows us to determine the dynamical regimes that 
can be expected in the HD 196885 system. 



cover the entire phase-space for this system, by exploring these 
two missing angles. 

Before studying the phase-space of the system, it is conve- 
nient to understand all the physical configurations. A scheme of 
the fundamental planes for the definition of the reference an- 
gles is shown in Figure [TJ For simplicity, we reserve for the 
inner orbit (planet) the index 1, and for the outer orbit (stel- 
lar companion) the index 2. The central star has mass mo=1.31 
M Q (Chauvin et al. 201 1), the planet has mass m\, and the stel- 
lar companion has mass ma, with semi-major axis a,, eccentric- 
ity e,, mean anomaly M,-, argument of pericenter w,, longitude 
of ascending node Q,, and inclination /,-. We also mark some 
additional angles that are useful for our study: the mutual in- 
clination i, the nodal longitude of the planet's orbit with re- 
spect to the binary's orbit Q, the difference of node longitudes 
AQ = Qi - Q2 in the plane of the sky, the argument of pericenter 
of the planet's orbit with respect to the binary's orbit w, and the 
angle Aw = u>\ — w. 

The dynamics of the system can be studied using different 
sets of the free parameters, the most important being (I\, AQ), 
(2, Q), and (2, w). The observed quantities obtained from fits to 
the data are deduced with respect to plane of the sky (I\ , Qi), but 
the dynamics of inclined systems is more adequately described 
using as reference the stellar companion (the orbit 2), since it 
is fully determined (Tab.[TJ and it is not much disturbed by the 
planet. 

All the reference angles are related between each other, so 
whatever is the choice that we adopt, one can easily determine 
the remaining angles by spherical trigonometry transformations 



3. 1 . Restricted quadrupolar problem 

We consider a binary star system composed of a primary (mo) 
and a secondary (2M2), and a massless planet (m\ =0) orbit- 
ing the primary star. The binary system's fixed orbit with pe- 
riod T%, semi-major axis «2 and eccentricity ez, is the natural 
choice of reference frame. The planet's osculating orbit with 
period T\, semi-major axis a\ and eccentricity e\, has orienta- 
tion with respect to the binary system's orbit defined by the an- 
gles i (relative inclination), w (argument of the pericen t re) an d 
Q (longitude of the ascending n ode). Following iKozail 0962), 
iKinoshita & Nakail d 19991 120071) . we write the Hamiltonian of 
the planet, expand up to quadrupole order in the semi-major axis 
ratio a.\lwi, and average with respect to the fast periods T\ and 
T2 , obtaining the secular quadrupole Hamiltonian 

F = C [(2 + 3 ef)(3 cos 2 i - 1) + 15 e\ sin 2 i cos(2w)] (7) 

witrfl 



C = 



Q 



m 2 



16(l- e 2)3/2 fl 3 



(8) 



The secular motion can be described by writing Hamilton's 
equations using Delaunay canonical variables 



(w,G = y^m fli (1 -<? 2 )] , 



(9) 



1 The expression for C (Eq. 7) in Kinoshita & Nakai (2007) should 
have a 1 , and not ai, in the nominator. 
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Fig. 2. Possible secular trajectories for the HD196885 system seen in the (to, i) plane (top), and in the (to, e\) plane (bottom). We 
show the trajectories in the restricted quadrupolar approximation (left) and using n-body numerical simulations of the real system 
(right). The dashed black curves in (a) separate zones of libration and circulation of the angle to. The black crosses in (a) and (c) 
indicate the location of Lidov-Kozai equilibrium points for this system. N-body numerical integrations in the (to, i) plane (b), and in 
the (to, e\) plane (d), preserve the same colors for the initial conditions in the restricted problem. 



Ul,H = Jgrnoai (1 -«?) cos?) , (10) 

where Q is the gravitational constant. Since the Hamiltonian 
(Eq.|7]) does not depend on Q, then the conjugate momentum H 
is constant. Moreover, as F (Eq. |7]i has only one degree of free- 
dom, to, we obtain secular trajectories by plotting level curves 
F - const, with 



H 2 



{&moa\) 



(1 



e,) cos i - const . 



(ID 



In our particular case (the HD196885 system), we know the 
current value of the planet's eccentricity (eio = 0.48) but we 
do not know the initial z'o or too- We can obtain possible secular 
trajectories for the planet by choosing values of it) and coq, then 
plotting the level curves of 



or 



F(i, to; h) = F(io, too; h) 



F(e\,to;h) - F(eid,(jL>(,\h) 



(12) 



(13) 



with h = (1 - e 2 ) cos iq. 

We show the secular trajectories described by expressions 
(fT2l and ( fT3l in Figure |2ja) and (c), respectively. Because of 



symmetry in the solution space we use the same color for tra- 
jectories with relative inclination i (prograde) or 180° - i (ret- 
rograde), and for librating trajectories around to = 90° or 
to - 270°. The separatrixes (dashed black curves in Fig.|2la)) 
mark the boundary between libration and circulation of the 
angle a> and are obtain ed by solving the implicit equation 
(iKinoshita&N akai 200^0 



F(i,to;h) = F s = 2C(3/i-l) 



(14) 



with h = (1 - e\ Q ) cos 2 /q for (a), i). 

In agreement with lLidovl(ll96lUl962l) and lKozail i 19621) . the 
regimes of secular motion consist of Lidov-Kozai equilibrium 
points (with i and e\ fixed) and Lidov-Kozai cycles which due 
to the conservation of h (Eq. fTTb exhibit coupled oscillations 
in i and e\. An orbit with e\ - 0.48 is at a Lidov-Kozai equi- 
librium point if (i = 47.2°, u = 90°), (i = 47.2°, o = 270°), 
(i = 132.8°, co = 90°) or (i = 132.8°, 10 = 90°). Lidov-Kozai cy- 
cles with co librating have h < 0.6 and F < F s (magenta, brown, 
cyan, yellow and pink orbits). Lidov-Kozai cycles with to cir- 
culating have h > 0.6 (red, blue and green orbits) or h < 0.6 
and F > F s (grey, orange and purple orbits). Moreover, high 
amplitude Lidov-Kozai cycles (grey, pink and orange orbits) can 



2 The solution to Eq. [14] is independent of h, as we would expect, 
since the separatrix is an invariant curve. 
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reach e\ near unity and could thus become unstable in the full 
problem. In particular, an orbit that passes through i - 90° must, 
by conservation of h (Eq. [TTJ, reach e\ — 1 (collision orbit). An 
example of this is the black orbit in Figure [2] which started at 
i = 90° and co = 90° (although it is not obvious from Fig.[2fa) 
due to Eq.Q~2]being undefined in this region). 

Note that Figure|2ja,c) are different from the standard Lidov- 
Kozai diagrams that show level curves of the Hamiltonian F at 
constant values h - Iiq. In our case eio is fixed but we show tra- 
jectories for different values of the initial inclination z'o, i.e., with 
different ho. In particular, this explains why librating orbits in the 
Lidov-Kozai regime do not encircle a Lidov-Kozai equilibrium 
point (magenta, brown, cyan, yellow and pink orbits in Fig. [2]). 
The location of the Lidov-Kozai equilibrium points in our di- 
agrams occur at the current planet's eccentricity (e\ - 0.48) 
hence circulating and librating orbits starting at co — 90° nor 
co - 270° will have eccentricity increasing (magenta, cyan, yel- 
low and pink orbits) or decreasing (red, blue, green and purple 
orbits) from this initial value (Fig. [2}. 



3.2. Numerical simulations 

In order to compare the orbital behavior in the full massive prob- 
lem with the restricted problem, we performed numerical inte- 
grations up to t - 30 kyr starting with the same initial con- 
ditions as in Figure Oa,c) (initial conditions from Table [TJ and 
co = 90° with i from 10 to 170 degrees, and co = 270° with i 
from 10 to 170 degrees). We show in Figure|2b,d) these numer- 
ical integrations as well as the position of the separatrixes and the 
Lidov-Kozai equilibrium points calculated with the quadrupolar 
approximation. We used the same colors as in Figure |2ja,c) to 
facilitate the comparison between the two models. 

We see that there is good agreement between the numerical 
integrations of the full problem (Fig.|3J),d) and the theoretical 
trajectories of the quadrupolar restricted Hamiltonian (Fig.|2^,c). 
The global behavior of the planet in the HD 196885 system 
is dominated by the dynamical regimes described above, i.e., 
Lidov-Kozai equilibrium points, librating and circulating orbits. 
In the full problem, the secular trajectories exhibit a drift from 
the secular solutions of the quadrupole Hamiltonian. This be- 
havior can be explained by i ncluding octupole and higher or- 
der terms in the Hamiltonian dFord et a l. 2000; Lask ar & Bouel 
2010; Lithwic k & Naozj[201 lb . The octupole Hamiltonian is not 
independent of Q and thus h (given by Eq. [TTJ is not a conserved 
quantity. Therefore, real secular trajectories (Fig. |2J,d) drift from 
the theoretical curves (Fig.|2^,c) due to changes in h. 

We also expect deviations from the restricted quadrupole 
model solutions if the planet's mass be comes comparab l e to th e 
star's masses although, according to Fara go & Laskarl (|2010), 
the topology of the phase space should not change. We see that 
this drift is more evident for orbits in the vicinity of the sep- 
aratrixes and orbits that reach e\ « 1 (nearly collision orbits). 
Note that the black trajectory in Figure |2lb,d) is the collision or- 
bit shown in Figure 0a,c) which is unstable after only 600 yr. 
Moreover, librating orbits around co - 90° seem to drift more 
than those around co - 270° which could be due to the high 
planet mass solutions in this region (Sect. [4] Fig. [3}. 

Also, the brown orbit in Figure |2jl is not a double loop struc- 
ture but two different initial conditions, one for (i < 90°) and the 
other for (i > 90°). In Figure[2j: both solutions coincide while in 
Figure [2ji those with i < 90° have a bigger loop than those with 
i>90°. 



4. Dynamical Analysis 

Several works on dynamics i n clos e bi nary systems exist. 
Among others iRabl & Dvorak! (1 19881) and iHolman & Wiegerj 
(1999), investigated the long-term stability of planets in copla- 
nar circular orbits near one of the stars. Expressions of critical 
semimajor axis for the planet (prograde orbits) were derived in 
function of the mass of binary components and eccentricity of 
the orbit. According to these approximations prograde orbits for 
this system are stable with semimajor axis less than ~ 3.8AU. 
Wiegert & Holman (1997) studied the stability of hypothetical 
terrestrial planets in the system cr-Centauri (whose semimajor 
axis and eccentricity are almost the same as in our system) for 
several mutual inclinations. However they fixed the longitude of 
the node and varied the semimajor axis (in our case the semi- 
major axis of the planet is very well established with the radial 
velocity technique). The authors also noticed that highly mutual 
inclined orbits are unstable. 

In this section we analyze the dynamics and the stability of 
the planetary system given in Table [TJ The best choice of vari- 
ables to study the global dynamics is the pair (/, co), since it al- 
lows to easily identify the fixed points and the different dynam- 
ical regimes (Sect.[3]l. However, this choice is not adequate as 
co\ imposes constraints on co, and the physical system has some 
restrictions in this frame (Eq.[2]or[5). Thus, we are left with the 
pair of angles (/, Q), which cover all the possible configurations, 
and are still independent of the observer. 

4.1. Orbital stability 

We constructed grids of initial conditions with 0.5° of resolu- 
tion and each point in the grid was then numerically integrated 
over 30 kyr using a Burlisch-Stoer based N-body code (precision 
better than 10~ 12 ) using astrocentric and osculating variables. 
During the integrations we computed the averaged MEGNO 
chaos indicator (Y): the regular orbits yield t o (Y) < 2, while 
larger values are indicative of chaotic motion (Cincotta & Simo 
2000). The MEGNO chaos maps uses a threshold that should 
be applied in order to avoid excluding stable orbits that did not 
converge to their theoreti cal value or tho se orbits that are weakly 
chaotic. Thus following Maffion e et al.l (1201 lh . the color scale 
shows "stable" orbits in blue up to (Y) ~ 2.5 (a particular choice 
based on integration of individual orbits for very long times and 
due to the characteristics of this system). 

Results are shown in Figure [3] We also show the grids in the 
observational frame (7;,AQ), and in the frame (i, co) for com- 
parison. We can see that highly mutually inclined systems (the 
strip within 70° < i < 1 10°) are highly unstable (this fact was 
already noticed in Wiegert & Holman ( 1997) for cr-Centauri but 
with much less resolution). 

In order to better understand the chaos regions we super- 
imposed contour level curves for mass of the inner planet (in 
Mjup)- As we change the values of (/, Q) we get corresponding 
value of I\ (Eq[T) and thus the real mass of the planet ranges 
between 2.98 Mj up (dashed lines) to 40 Mj llp (with larger val- 
ues confined to little circles). Inside these circles the mass of the 
planet tends to infinity. Mass level curves become straight lines 
in the observational frame, since m\ <x 1 / sin I\, but the dynami- 
cal regimes are not perceptible, except instability corresponding 
to high masses for I\ - 0° and 180°. 

In the frame (/, co), all the main structures become under- 
standable: the system is unstable for high values of the mutual 
inclination, for high masses, but also around the separatrixes be- 
tween dynamical regimes (Sect. [3). In addition, we observe that 
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some initial conditions are not possible near high mutual incli- 
nations (i ~ 90°). The reason for this "forbidden" zone is that 
some initial conditions in the frame (i, Q) are degenerate and 
correspond to identical solutions in the frame (i, a>), although 
with different values for the mass of the planet. 

The size of the "forbidden" area depends on the inclination 
I2 of the stellar companion to the plane of the sky (Eqs.[TJ[2]i. For 
h - 90° there is a perfect correspondence between the frames 
(i, Q) and (i,u>), while for I 2 - 0° or 180°, there is a global 
degeneracy. For the HD 196885 system, we have I2 - 116.8°, 
that is, about 27° above 90°. As a consequence, the size of the 
degenerate area is about +27° around i - 90° in the frame (i, to). 
The delimitation of this "forbidden" zone also coincides with 
the areas of large masses (Fig.[3];,d). In the particular case of the 
HD 196885 system, in the areas with overlapping solutions, one 
always corresponds to unstable orbits (the one that is close to the 
borders of the zone). Therefore, for simplicity, we can merge the 
two figures in a single diagram, ignoring the areas with i < 90° 
in Figure [3jc), and i > 90° in Figure [3jd)- F° r the remaining of 
the paper we always adopt this procedure in the frame (i, a>). 

We used the MEGNO chaos indicator because it allows us to 
rapidly distinguish between regular and chaotic orbits. MEGNO 
is very useful in these type of studies because we do not need 
to check every initial condition for long periods of time. We 
can just follow a set of orbits in representative regions of the 
MEGNO map to understand how chaos leads to instability. 
Although the MEGNO indicator only provides indications on the 
regularity of the orbits, we have verified that the chaotic orbits 
marked in red are actually unstable in the sense that the planet 
reaches distances to the central star smaller than the radius of the 
star itself (we set the limit at 0.005 AU which is the radius of the 
Sun, an overestimate quantity because for this system, the cen- 
tral star has 1.3 ^sun)- The timescales of this instability depend 
of the chaotic region that is considered: 

1. Around the Kozai separatrices (red regions), chaotic orbits 
become unstable in short timescales (less than ~ 25000 
years). 

2. When the initial mutual inclination is between 80° < i < 
100° the orbits are very unstable and the planet collides with 
the central star in less than 1000 years. As we explained in 
Section[3] an orbit that passes through i - 90° must, by con- 
servation of h, reach e\ - 1 (collision orbit). 

3. Near the edge of the forbidden region: (i ~ 60°, u> < 180°) 
and (/ ~ 118°,w < 180°), the planets have masses > 
AQMj llp {i\ ~ 0°), thus the strong interactions with both stars 
make the planet either collide with the central star or escape 
from the system. 

4. There is a wide region of unstable orbits with high values of 
MEGNO inside the left-hand Kozai separatrix (55 < i < 60° 
and 45° < u> < 125°) where the collision times goes from 
10 5 to 10 7 years. 

Within the Kozai libration islands around to - 270° the 
stable initial conditions have eccentricities that remain always 
lower than 0.90. Inside the Kozai libration islands around a> - 
90°, the retrograde conditions show the same behaviour, while 
some prograde conditions only survive for times from 10 7 to 10 8 
years. Outside the separatrix, in the prograde region, the chaotic 
orbits (in green at bottom left hand of the graph, where (Y) > 2.5 
and 4° < i < 10°) are stable not reaching ever eccentricities 
above 0.7. 

Long-term numerical simulations showed that the regions 
marked in blue/green in the MEGNO map do not present de- 



tectable dynamical instability, at least in time-scales of the or- 
der of 10 9 years. Due to the high mass ratio, the effect of high 
order mean-motion resonances (MMR, hereafter) may be non- 
negligible (the period ratio is 20.95). In fact, there are a multi- 
tude of high order MMR whose overlap could be responsible for 
the slow chaos regions. However, these high order MMR are dif- 
ficult to identify, and their relative importance will depend on the 
exact orbital elements of the system. Even secular effects could 
be the explanation for the slow chaos regions in green. In partic- 
ular, we refer that chaos nearby the Kozai separatrix has already 
been observed in the octupole problem (Lithwick & Naoz 201 1 ). 
We prefer to conclude that the best representation of the sys- 
tem corresponds to a region of regular motion (the region marked 
in blue). Although this may seem arbitrary, it is important to re- 
call that so far, in any known planetary system, there are no giant 
planets displaying significant chaotic motion. Secondly, regions 
of regular motion are expected to be more robust with respect to 
additional perturbations and planetary formation. 

4.2. Dynamical regimes and regular orbits 

MEGNO is a very efficient tool to identify chaotic motion, but 
it is not suited to distinguish between different types of regular 
orbits, that is, orbits that are more or less unperturbed, or orbits 
that experience intense secular variations. In order to evaluate 
the orbital behavior, in Figure |4] we show the amplitude of the 
inclination variations for the three frames, together with the sep- 
aratrixes between dynamical regimes obtained in Sect. [3] 

We observe that unstable regions (Fig. [3} mainly coincide 
with the regions of large amplitude variations of the mutual in- 
clination. Indeed, according to expression (fill and as discussed 
in Sect. [3] large variations in the inclination up to i « 90° can 
increase the eccentricity of the inner planet to very high values, 
which gives rise to close encounters with the stellar companion. 

The more regular orbits of the system are those which are 
close to the Lidov-Kozai equilibrium points, and those corre- 
sponding to coplanar orbits. Also notorious is that retrograde 
orbits (i > 140°) are in general less perturbed, probably be- 
cause close encounters last less time. As we have seen before 
(Fig. [3), the zones of higher mass for the planet introduce some 
instability, but the global picture is dominated by the dynamical 
regimes: trajectories that are close to the separatrix and nearly 
collision orbits are clearly the most unstable (cf. Sect. 3). 

4.3. Eccentricity of the inner planet 

The eccentricity of the inner planet is a key variable to under- 
stand the stability and the evolution of the system, since high val- 
ues lead to close to collision with the central star. Indeed, while 
the outer companion may destabilize the orbit of the planet, en- 
counters with the central star may give rise to tidal effects and a 
subsequent secular evolution of the orbit (Sect. l4.5l l. In addition, 
among all orbital parameters listed on Table [TJ the eccentricity 
<?i is the only one for which secular modifications can be ob- 
served using radial velocity data (the argument of the periastron 
cj\ also varies, but it is not directly related to the physical or- 
bit). In Figure |5]we show the secular period associated with the 
largest eccentricity oscillation of the inner planet (in years), the 
amplitude of these eccentricity oscillations, and the maximum 
eccentricity, e\, attained during the integration. 

As discussed previously, while the behavior of high inclina- 
tion orbits (i.e. Lidov-Kozai cycles) is well described with the 
quadrupole Hamiltonian, the secular evolution of nearly copla- 
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Fig. 3. MEGNO chaos indicators for different sets of variables, (a) (i,Ci), (b) (/i, AQ), (c) (/, w) with 0° < Q < 180°, and (d) (i,cj) 
with 180° < Q. < 360°. Labelled lines give the mass of the planet in Mj up . Orbits can be considered stable for (Y) < 2. 



nar orbits must take into account the octupole Hamiltonian. In 
particular, the eccentrity of co planar orbits chang es secularly 
solely due to the octupole term (Lee & Peale 2003). The eccen- 
trity of low to medium inclination orbits is modulated by the oc- 
tupol e and quadrupole secular periods (Krymolowski & Mazeh 
1999). Moreover, the quadrupole secular period is typically 
shorter than the octupole secular period (Fig. |6). Octupole and 
higher order terms are also responsible for making the orbit s near 
the separatrixes chaotic Lithwick & Naozl (1201 II) . 

In Fig. |5|a) we see the secular period associated with the 
largest eccentricty, e\, variations. From Fig. Ha) we see that 
retrograde orbits have secular periods longer than direct or- 
bits. Among stable areas, medium amplitude Lidov-Kozai cycles 
have the largest eccentricity variations (Fig.[6]and Fig.|3J;,d) and 
periods of around 1000 yr (Fig. |6]and Fig. 0a). The eccentricity 
oscillations of small amplitude Lidov-Kozai cycles occur on a 
1500 yr timescale while nearly coplanar orbits also exhibit small 
oscillations in e\ but on a timescale of around 5500 yr (Fig. |6j. 
Although we may be able to measure eccentricity variations in 
the radial velocity data (for instance, by detecting a drift in the 
data), in practice, we will have only a short observation timespan 
(typically a few years). Therefore, it will be difficult to distin- 
guish small amplitude Lidov-Kozai cycles from nearly coplanar 
orbits solely based on the different secular timescales (e.g. com- 
pare magenta and red orbits in Fig. |6). We may, however, be able 
to identify if the planet is on a medium amplitude Lidov-Kozai 



cycle due to the large eccentrity variations on a short (1000 yr) 
timescale. 

The amplitude of the eccentricity oscillations (Fig. [2) does 
not differ much from the amplitude of the mutual inclination os- 
cillations (Fig.|4j;). This behavior was expectable since these two 
variables are correlated (Eq.[TTV As for the inclination, we can 
easily identify the different dynamical regimes. 

In Figure [5};, we plot the quantity log I0 (l. - ei,„ MA ), where 
e\,max is the maximal excentricity attained by the planet during 
the integration. Thus a value of log 10 (l. - e\ jnax ) - -3 (in red) 
represents an eccentricity of e\ jm]x - 0.999. Light-blue regions are 
those with e\^ nax ~ 0.9, dark-blue for those with ei. mflA ~ 0.7, and 
finally light-gray for those where e\ max ~ 0.5. The initial condi- 
tions within violet region inside the Kozai-Lidov separatrix can 
survive for times from 10 8 to 10 9 years. Note that this confirms 
that the orbits marked in red in MEGNO map collide with the 
central star because e ~ 1 (as described in Subsection l4.U 

4.4. Different initial conditions 

Since the orbit al parameters i n TableQ~]contain some undetermi- 
nations (IChauvin et al.ll201ll) . we also explored the stability of 
the system within one cr of the best fit parameters. In particular, 
we wanted to test the impact of variations in the outer compan- 
ion, which is the most unconstrained. For that purpose, we have 
made grids at the upper limits at - 21.86 AU and ei - 0.45, and 
at bottom ones a^ - 20.14 AU and e2 = 0.39 as in Figure|3] We 
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Fig. 4. Amplitude of the mutual inclination for different sets of 
variables, (a) (/, Q), (b) (I\, AO), (c) (i, oj). Centers of the Lidov- 
Kozai equilibria are marked as black circles and the separatrix is 
calculated with the quadrupolar approximation (Sectj3)- 
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do not show the results, since unstable (red regions) are located 
in same positions (i ~ 90°, and around the Kozai separatrix). 
Only regions of moderate chaos (green regions) varied their po- 
sitions and size a little bit. However, those regions are still pos- 
sible solutions. Therefore, we conclude that the global picture in 
Figure|3] does not change much for different sets of initial con- 
ditions, that is, all the conclusions in this paper are still valid 
within the errorbars of the published parameters. 



We also make a grid for the restricted case (i.e. the mass of 
the planet m\ - 0) using the best- fit parameters from TableQ] 
All the chaotic green regions disappeared and we were left only 
with regular regions (blue) and close encounter regions (in red). 

4.5. Tidal evolution 

For some orbital configurations, the eccentricity of the inner 
planet may reach very high values, and become close enough 



C. A. Giuppone et al.: Dynamical analysis and constraints for the HD 196885 system 




2000 4000 6000 8000 10000 
Time [years] 




2000 4000 6000 8000 10000 
Time [years] 



Fig. 6. Time variation of eccentricity for the orbits in Fig. |2] 
showing the two secular periods. The color code is the same than 
Fig. [2] Initial conditions outside the separatrix (top), and inside 
the separatrix (bottom). 



to the central star at periastron to undergo tidal effects. For an 
unperturbed orbit, the secular evolution of the eccentricity by 
tidal effect can be given by dCorreiall2009l) : 



e\ = -Kf(e\)e\ 



with 



/(«) 



1 + 
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(15) 



(16) 



(17) 



where ki is the second Love number, Q is the tidal dissipation 
factor, R is the radius of the planet, and «o = Oi(l - e\) - Const. 
The solution of the above equation is given by 
(Correia&Laskar2010) 



F(ei) = F (e initial) exp(-Ki) . 



(18) 



where F(e) is an implicit function of <?, which converges to zero 
as t — * +03. The characteristic time-scale for fully dampening 
the eccentricity of the orbit is then r ~ 1 IK. The above expres- 
sion is only valid for unperturbed orbits, but in the HD 196885 




l, max 



Fig. 7. Damping time for the eccentricity of the inner planet 
(with t = l/K, Eq.EJ. For Q ~ 1000, only for e hmax > 0.96 
the eccentricity is damped within 2 Gyr, the estimated age of the 
system (j/Q = 2 x 10~ 3 Gyr, red line). 



system the initial eccentricity can be seen as the maximal eccen- 
tricity attained on a Kozai cycle, e\ ttnax , since it is at this point 
that dissipation is maximized. Therefore, r provides a minimal 
estimation of the dampening time. For the orbits for which the 
eccentricity can be damped during the age of the system, we can 
exclude them from the possible orbital parameters of the planet, 
since the present eccentricity is still near 0.5 (Tab.Q}. 

In Figure [7] we show the evolution time-scale for differ- 
ent initial maximal eccentricities, adopting k.2 = 0.5 and R = 
1 .2,Rj up . We observe that for Q ~ 1000, a value similar to Jupiter, 
only for ei. mflA > 0.96 the eccentr icity is damp ed within 2 Gyr, 
the estimated age of the system dCorreia et al. 2008). Even for 
Earth-like planets (Q ~ 10), only e^ max > 0.92 can be dissipated 
during the same amount of time. Since e\ > 0.96 corresponds to 
unstable orbits, we conclude that tidal effects cannot be used to 
exclude any of the stable initial conditions compatible with the 
observational data (Tab.[T]i. 



5. Formation hypothesis 

From the formation point of view, the orbital excitation ex- 
erted by Lidov-Kozai cycles implies destructive, high-velocity 
collisions among planetesimals, inhibiting formation of massive 
objects (lLissauerill993h . However the Lidov-Kozai equilibrium 
gives itself a protection mechanism to get stable configurations 
when mutually inclined bodies are considered. 

Thebault (201 1) published a deep study of the possible sce- 
narios of formation for the planet in the HD 196885 binary 
system. The author considered axisymmetric static gas disc 
(no self gravity) and estimated the impact velocities amongst 
a population of circumprimary planetesimals. A main conclu- 
sion was that the circumprimary disc is strongly hostile to plan- 
etesimal accretion, especially the region around 2.6 AU (the 
planet's present location) where the binary perturbations induce 
planetesimal-shattering velocities of more than 1 km/s. The re- 
gion around 2.6 AU is strongly hostile to planetesimal accretion, 
even for highly inclined orbits and alternative solutions were 
proposed to justify their existence ( e.g., different initial co nfigu- 
ration of binary or disk instability). lOuintana et al.l (12002) stud- 
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ied the formation in a - Centauri system using a disk of com- 
bined large and small bodies. They concluded that it is possible 
to form planets when an inclined disk is considered, although 
~ 95% of the initial mass is lost by the end of the simulation. 
Their results predict a wider region to form more retrograde or- 
bits (/ = 180°) than prograde ones and that maybe large plane- 
tary embryos could be form near the star (from 0.5 to 1.5 AU) 
in inclined orbits. However as the initial semimajor axis distri- 
bution has an upper limit set at ~2 AU, this may bias the final 
results for the allowed regions. 

However, Batygin e t alj d2011l) address the formation for 
wide binary systems (abinary ~ 1000 AU), concluding that is pos- 
sible to form a single planet in an inclined orbit (Lidov-Kozai 
regime), if taking into account the self gravity of proto-planetary 
disk (that means planetesimals embedded in the gaseous disk). 
They also pointed out that the evolutionary process of forma- 
tion of a planetary system at the Lidov-Kozai equilibrium is 
necessarily non-unique and proposed a scenario were a multi- 
ple system could be formed protected by Lidov-Kozai cycles 
and subsequent instabilities remove all the remaining planets. 
Unfortunately no mention to close binary systems with this ap- 
proach could be found. 

The secular dynamics of small planetesimals play a funda- 
mental role in establishing the possibility of accretional colli- 
sions in such extreme cases. The most important secular pa- 
rameters are the forced eccentricity and the secular frequency, 
which depend on the initial conditions of the particles, as well as 
on the mass and orbital parameters of the sec ondary star (e.g. 
Thebault 201lHBeauge et alJl20ld) . However iGiuppone et alJ 
( 20 1 1 ) pointed out that for these kind of compact systems, some- 
times the frequencies are not quite well determined from first or- 
der approximations (second order on mass are needed), and thus 
probably, some works on formation based in these approxima- 
tions should be reviewed. 

As for the giant planets in the Solar System, that are sup- 
posed to have migr ated owing to their i nteraction with a disk of 
planetesimals (e.g. lTsiganis et all2 005), we may assume that the 
same occured with the orbit of the planet in the HD 196885 sys- 
tem. Using the secular model described in ICorreia et al.l (1201 1) 
for the full system (two stars and one planet), we have performed 
some simple experiments on the early evolution of the system, 
by migrating hypothetical initial orbital parameters of the planet 
into the present ones. 

We have independently tested the evolution of the semi- 
major axis and t he eccentricity of the inner orbit, u sing an ex- 
ponential decay dBeauge e t al. 20061 iLee et"alll2007h : 



«i = — -exp(-?/-r m ) 



e\ 



T„,/10 



(19) 



with t,„ = lOMyr. For the semi-major axis, there is no varia- 
tion in the orbital configuration, since the ratio 01/02 is a fac- 
tor in the quadrupolar Hamiltonian (Eq.|8), so it only modifies 
the time-scale of the evolution. However, for the eccentricity 
we observe that all trajectories which begin inside the libra- 
tion zone migrate into one of the Lidov-Kozai equilibrium po- 
sitions (Fig. [8}. Trajectories in the circulation zone remain more 
or less unchanged, only the amplitude of the mutual inclination 
is damped. Therefore, we can assume that if the planet in the 
HD 196885 was able to form in libration, it will be presently ob- 
served near the Lidov-Kozai equilibria. 
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Fig. 8. Early evolution of the planet, when the eccentricity is 
damped from higher values due to the interaction with a disk 
of planetesimals. For each simulation, the color of the position 
of the planet becomes darker with time. In the circulation zone 
only the amplitude of the mutual inclination is damped, but for 
orbits starting in the libration area, the planet evolves into the 
Lidov-Kozai equilibria. 



6. Discussion and conclusions 

We have studied the dynamics of the HD 196885 system using 
the present determination of the orbital parameters (Tab.[T} as 
starting point. We developed an insight full analysis of the 3-D 
space exploring the free parameters and choosing different sets 
of angular variables. We found that the most likely configura- 
tions for the planet in the HD 196885 system is either nearly 
coplanar orbits (prograde and retrograde), or highly inclined or- 
bits near the Lidov-Kozai equilibrium points, i - 90° = +47°. 

Among coplanar orbits, the retrograde ones appear to be less 
chaotic, possibly because close encounters last less time than 
in the prograde case. For the orbits at the Lidov-Kozai equilib- 
rium, those around u> - 270° are more reliable than those around 
90°, since the mass of the planet is smaller in the first situation. 
Present formation scenarios for this kind of systems are unable 
to distinguish between the two possibilities (coplanar or Lidov- 
Kozai equilibrium). However, a simple simulation on the evolu- 
tion of the initial eccentricity of the inner orbit shows that if the 
planet was able to form in the libration zone, it will evolve into 
the Lidov-Kozai equilibrium point. We also tested the effect of 
tides, but they are too weak in the HD 196885 system, although 
they should be taken into account for tighter systems of this kind. 

Although there is a wide variety of stable initial conditions 
for the planet in the natural frame defined by the orbit of the two 
stars (/, Q), from the observer's point of view (plane of the sky) 
there are some restrictions. Indeed, looking at Figure 0b) we see 
that stable areas occur around: 






Q! ~ 80° 
65° kozai prograde 

115° coplanar prograde 



Qj ~ 260° 

coplanar retrograde 

kozai retrograde 



Thus, by continuing to observe the system one will be first able 
to determined the dynamical regime of the planet (coplanar or 
Lidov-Kozai equilibrium) and later its exact position. Moreover, 
it may be possible to identify medium amplitude Kozai cycles 
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due to their large eccentricity oscillations on 1000 yr timescale. 
Small amplitude Kozai cycles and nearly coplanar orbits have 
small eccentricity oscillations on distinct timescales (1500 yr 
and 4000 yr timescales, respectively). However, in practice, it 
may be difficult to distinguish these two configurations as we 
are limited to short observation timespans. 
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